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ABSTRACT 

This  thesis  estimates  the  frequency  response  of  a  network  where  the  only  data  is  the 
output  obtained  from  an  Autoregressive-moving  average  (ARMA)  model  driven  by  a 
random  input. 

Models  of  random  processes  and  existing  methods  for  solving  ARMA  models  are 
examined.  The  estimation  is  performed  iteratively  by  using  the  Yule-Walker  Equations 
in  three  different  methods  for  the  AR  part  and  the  Cholesky  factorization  for  the  MA 
part.  The  AR  parameters  are  estimated  initially,  then  MA  parameters  are  estimated 
assuming  that  the  AR  parameters  have  been  compensated  for.  After  the  estimation  of 
each  parameter  set,  the  original  time  series  is  filtered  via  the  inverse  of  the  last  estimate 
of  the  transfer  function  of  an  AR  model  or  MA  model,  allowing  better  and  better  esti- 
mation of  each  model's  coefficients.  The  iteration  refers  to  the  procedure  of  removing 
the  MA  or  AR  part  from  the  random  process  in  an  alternating  fashion  allowing  the 
creation  of  an  almost  pure  AR  or  MA  process,  respectively.  As  the  iteration  continues 
the  estimates  are  improving.  When  the  iteration  reaches  a  point  where  the  coefficients 
converse  the  last  VIA  and  AR  model  coefficients  are  retained  as  final  estimates. 
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I.     INTRODUCTION 

In  this  thesis,  an  iterative  approach  is  presented  to  estimate  the  frequency  response 
of  a  network  where  the  only  data  is  the  output  obtained  from  an  Autoregressive-moving 
average  (ARMA)  model  driven  by  a  random  input.  The  AR  parameters  are  estimated 
initially,  then  the  MA  parameters  are  estimated  assuming  that  the  AR  parameters  have 
been  compensated  for.. 

To  find  the  frequency  response  of  an  ARMA  network  two  problems  have  to  be  ad- 
dressed. One  is  due  to  the  shortness  of  the  observed  data  (time  limitation)  and  the  re- 
sulting distortion  of  the  estimated  correlation  function.  The  second  problem  is  due  to 
the  nonlinear  combination  of  the  coefficients  of  the  MA  part  as  they  appear  in  the  esti- 
mate of  the  correlation  function. 

To  achieve  this,  we  estimate  coefficients  of  MA  part  and  AR  part  iteratively.  We 
use  the  Yule-Walker  Method  for  AR  coefficients  estimation  and  the  Cholesky  decom- 
position for  MA  part  estimation.  We  assume  that  the  ARMA  network  output,  i.e.,  the 
observed  signal  is  produced  by  white  Gaussian  noise  driving  the  network.  Three  differ- 
ent methods  are  used  to  estimate  the  AR  coefficients.  To  minimize  the  effect  of  the 
MA  part  in  the  correlation  function,  correlation  lags  greater  than  the  correlation  length 
of  the  MA  part  are  used  in  estimating  the  initial  AR  coefficients.  The  iteration  refers 
to  the  procedure,  which  removes  estimated  MA  or  estimated  AR  contribution  from  the 
output  of  the  ARMA  model  in  an  alternating  fashion.  This  allows  better  and  better  es- 
timates of  the  AR  and  MA  coefficients  as  the  iteration  continues.  For  the  iterative  AR 
coefficient  estimation  three  methods  can  be  used,  denoted  by  method  1,  method  2  and 
method  3.  Method  1,  2  and  3  use  correlation  lags  starting  at  p+  1,  0,  and  0  respectively, 
where  the  first  two  methods  use  a  square  matrix  inverse  and  the  third  method  uses  a 
Pseudo  matrix  inverse. 

Introduction  to  models  of  random  processes  is  presented  in  Chapter  II.  Existing 
techniques  to  solve  for  the  parameters  of  ARMA  models  are  presented  in  Chapter  II. 
Chapter  IV  includes  simulation  results.  Simulation  studies  employed  the  Matlab  pack- 
age on  the  IBM  PC/AT.  Conclusions  and  recommendations  are  presented  in  Chapter 
V. 


II.      MODLLS  OF  RANDOM  PROCESSES 

A.      AR  PROCESS  MODEL 

We  say  that  x(n)  is  an  autoregiessive  process  of  order  p,  or  simply  an  AR(p)  process 
if  it  satisfies  the  difference  equation, 


x(n)  +  ci\x(n  —  J  j  +  ....  +  (ipx(n  —  p)  =  i(n) 


(2.1) 


where  alt  a2.  ••■-,  «r  arc  constant  cocilicients.  and  e(n)  is  a  pure  random  process  [Rcf.  1]. 
Equation  (2.1)  may  be  written  in  the  following  form, 


*{»)  =  ~  i_,aiA"  ~  k)  +  E(") 
A  realization  of  L:q.(2.2)  is  illustrated  in  Figure   1. 


(2.2) 


Figure   1.       Autoregiessive  Model  of  Order  p 


The  system  transfer  function  II(z)  between  the  input  c.(n)  and  the  output  x(n)  for 
the  AR  model  shown  in  Figure  1  is 


H{z)  =  — 7—  = ; l- (2.3) 

Az)        l+a]z-]+a2z-2  +  ....  +  apz-p 

It  is  required  that  A(z)  has  all  its  roots  within  the  unit  circle  in  the  z-plane  which 
guarantees  that  H(z)  is  a  stable  and  causal  filter. 

The  form  of  equation  (2.3)  illustrates  that  AR  models  have  finite  poles  but  no  ze- 
roes. Hence,  this  model  is  sometimes  called  an  All-Pole  model.  Because  an  AR  model 
has  no  poles  outside  or  on  the  unit  circle,  it  has  the  strict  minimum-delay  property  and 
hence  is  always  invertible.  In  general,  minimum  delay  means  that  the  transfer  function 
must  have  no  poles  outside  the  unit  circle,  but  can  have  poles  on  the  unit  circle.  Strictly 
minimum  delay  means  that  the  transfer  function  has  no  poles  outside  or  on  the  unit 
circle. 

The  AR  model  is  also  called  an  Infinite  Impulse  Response  (II R)  filter.  According 
to  definition  (2.2),  output  x(n)  depends  on  past  values  of  the  output  and  on  the  present 
input.  Because  of  this,  it  is  also  referred  to  as  a  Pure  Feedback   system. 

The  power  spectral  density  for  AR  models  is  given  by 

'^■TwF  {2A) 

where 


A(j)=\+ Yfk*~jnfkT 


and  T  is  the  sampling  interval  [Ref.  2]. 

B.      MA  PROCESS  MODEL 

The  sequence  x(n)  is  said  to  be  a  moving  average  process  of  order  q  (denoted  by 
MA(q))  if  it  satisfies  the  difference  equation, 

x{n)  =  b0s{n)  +  b^{n  -  1)  +  ....  +  bqz{n  -  q)  (2.5) 

where  b0,  bu  ....,  6?are  coefficients,  and  e(n)  is  a  pure  random  process  [Ref.  1]. 
Equivalently,  we  may  write, 


/<=0 


(2.6) 


We  may  say  that  the  output  of  an  MA  model  depends  only  on  present  and  past 
values  of  the  input,  i.e.,  there  is  no  feedback  in  an  MA  model. 
A  realization  of  Eq.  (2.6)  is  illustrated  in  Figure  2. 
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Figure  2.       Moving  Average  Model  of  Order  q 


The  system  transfer  function  H(z)  between  the  input  z(n)  and  the  output  x(n)  for 
the  MA  process  is 


H{z)  =  B(z)  =  bQ  +  bxz  l  +  b2z     +  ....  +  b,i 


(2.7) 


This  transfer  function  has  q  finite  zeroes,  but  no  poles.  Hence,  the  MA  model  is 
also  called  an  All-Zero  model.  MA  models  are  invertible  if  and  only  if  B(z)  has  no  zeroes 
outside  the  unit  circle,  nor  on  the  unit  circle. 

MA  models  are  also  called  Finite  Impulse  Response  (FR)  filters. 

The  power  spectral  density  for  an  MA  model  is  given  by 


PuA(f)  =  o2t\B(f)\2 


(2.8) 


where 


B(J)  =  b0  +  YJhe~j7:fkT 


k=\ 


and  T  is  the  sampling  interval. 

C.      ARMA  PROCESS  MODEL 

We  say  that  x(n)  is  an  autoregressive-moving  average  process  of  order  (p,q)  or 
simply  an  ARMA(p,q)  process  if  it  is  satisfies  the  difference  equation, 

x(n)  +  a}x(n  —  1)  +  ....  +  OpX(n  —  p)  =  bQe{n)  +  bxz(n  —  1)  +  ....  +  bqe(n  —  p)       (2.9) 

where,  again,  au  ....,  ap,  b0, ....,  bq  are  coefficients  and  z[n)  is  a  pure  random  process  [Ref. 

!]• 

Equation  (2.9)  may  be  written  as, 

P  <7  oo 

x(„)  =  -  ^a^in  -  k)  +  ^Jbkz{n  -  k)  =  ^\kz{n  -  k)  (2.10) 

k=\  k=0  A-=0 

The  assumption  b0  =  1  can  be  made  without  any  loss  of  generality  because  the  input 
c(n)  can  always  be  scaled  to  account  for  any  filter  gain  [Ref.  2]. 
A  realization  of  Eq.  (2.10)  is  illustrated  in  Figure  3. 
The  system  transfer  function  for  the  ARMA  process  is  given  by 

B{z)        1  +  b,z~l  +  b2z~2  +  ....  +  baz~q 

H(z)  =  -j-{  = ~ q—  (2.11) 

A(z)        1  +  axz    X  +  a2z  2  +  ....  +  apz  p 

where  A(z)  is  the  z-transform  of  the  AR  part  and  B(z)  is  the  z-transform  of  the  MA  part. 

Both  polynomials  A(z)  and  B(z)  are  assumed  to  have  all  of  their  zeros  within  the 
unit  circle  of  the  z-plane  to  guarantee  that  H(z)  is  a  stable  minimum-phase  invertible 
filter. 

The  power  spectral  density  for  ARMA  models  is  given  by  [Ref.  2] 
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Figure  3.       Autoregressive-Moving  Average  Model  of  Order  (p,q) 


D.      RELATIONSHIPS  OF  RANDOM  PROCESSES 

The  Wold  decomposition  theorem  [Ref.  3]  relates  the  AR,  MA  and  ARM  A  models. 
It  shows  that,  if  the  Power  Spectral  Density  is  purely  continuous,  any  AR  or  ARMA 
process  can  be  represented  by  a  unique  MA  model  of  infinite  order. 

Another  important  theorem  which  is  stated  by  Kolmogorov  [Ref.  4]  says  that  any 
ARMA  or  MA  process  can  be  represented  by  an  AR  process  of  infinite  order. 


To  illustrate  these  theorems,  we  model  [Ref.  5]  an  ARM A(  1,1)  process  by  an 
AR(oo)  or  by  an  MA(oo)  process.  From  equation  (2.1 1),  the  system  transfer  function  for 
the  ARM A(  1,1)  process  is 


H(z)  = 


1+ V~1 
1  4-  a,z 


If  we  use  AR(oo)  process  to  represent  ARMA(1,1)  process,  where 

H&  = -l     1       -2 

1  +  CiZ        +  C-yZ        +  .  .  . 


and 


00  -i 

-k       1  +  a\z 


c(Z) = i + Y,^~k = - 


+  b,z~] 


k=\ 

By  using  synthetic  division  we  find  that 

C(z)  =  1  +  fa  -  bx)z^  +  {b]  -  axbx)z~2  +  {b\  -  axbx)z~3  +  .  .  . 
Hence  inverse  z-transform  of  az~m  is  ad(k  —  m)   [Ref.  6]  and 

fl        if    k  =  0 

(o       else, 

the  inverse  z-transform  of  C(z)  is 

ck  =  S(k)  +  (ax  -  bx)5(k  -  1)  +  {b2x  -  axbx)d(k  -  2)  +  [b]  -  axb2x)S(k  -  3)  +  .  .  . 


or 


1  if  k  =  0 


Cfr  = 


k      l^-b^-b/-1  if/c>l 

If  we  use  a  finite  order  AR(p)  we  should  choose  p  to  satisfy  c,+1~0  or,  equivalently 

Therefore,  a  high-order  AR  model  will  be  required  when  the  zero  of  the  ARMA 
process  gets  closer  to  the  unit  circle. 


In  a  similar  way  if  we  use  an  MA(oo)  process  to  represent  an  ARM  A(  1,1)  process, 
let 

H{z)  =  d0  +  dxz~x  +  d2z~2  +  ... 

where 

00  -1 

V  j    -k       !  +  V 


D{z)  =  2J#     = 


1  +  a,z 

k=0  ' 


By  using  synthetic  division  as  we  did  above,  the  inverse  z-transform  of  D(z)  will  be 

j  1  if  k  =  0 

*"!(*,-«,)( -a,)*"1      if  k>\ 

If  we  use  a  finite  order  MA(q)  model,  we  should  choose  q  to  satisfy  dq_^0  or, 
equivalently  ofcO. 

Therefore,  a  high-order  MA  model  will  be  required  when  the  pole  of  the  ARMA 
process  gets  closer  to  the  unit  circle. 

E.      RELATIONSHIP  OF  AR,  MA,  AND  ARMA  PARAMETERS  TO  THE 
AUTOCORRELATION  SEQUENCE 

In  this  section,  we  will  present  the  relationship  of  the  model  parameters  to  the 
autocorrelation  sequence  [Ref  2]. 

If  we  multiply  Eq.  (2.10)  by  x*(n  —  m)  and  take  expectation,  the  result  will  be 

p  1 

E{x(n)x*(n  —  m)}  =  —   /  akE{x(n  —  k)x*(n  —  m)}  +  /  pkE{z(n  —  k)x*(n  —  m)}  (2.13) 

fc=1  k=0 

where  the  superscript  *  is  used  to  denote  the  complex  conjugation. 
Equation  (2.13)  may  be  written  as, 

p  <7 

rxx(m)  =  ~  2_jakrxxim  ~k)  +  ZJlfJ™  ~  k)  (2- 14) 

k=\  k=0 

The  cross  correlation  between  the  input  and  the  output  can  be  written  as, 


S 


rjl)  =  E{z{n  +  f)x*(n)}  =  E{  e(n  +  I) 


B*(#l)  +  2//«*(»-*) 


k=\ 


«M-'JO  +  XVr«('+*) 


(2.15) 


k=\ 


where  A0  =  1  by  definition  (Eq.  2.10). 

If  we  assume  that  the  driving  sequence  is  a  white  noise  process  of  zero  mean  and 
variance  a]  then  [Ref.  2] 


rM  = 


0 


for    />0 
for    /  =  0 


(2.16) 


efji*^  for    /<0 


When  we  substitute  Eq.  (2.16)  in  Eq.  (2.14),  we  get  final  relationship  between  the 
ARMA  parameters  and  the  autocorrelation  sequence. 


rxx(m)  = 


rxx*(  —m)  for    m<  0 

P  9 


~  2]<ikrxx(m  ~  k)  +  °]  Yu  bkh*\k-m)  for     0<m<q 

k=\  „  k=m 

-  Z^ak>'xx(m  ~  k)  for    m  >  q 


k=\ 


The  relationship  between  the  autocorrelation  sequence  and  a  pure  autoregressive 
model  may  be  written  by  setting  q  =  0  in  Eq.  (2.17) 


rxx(m)  = 


~  2^akrxx(m  -  k)  for      m  > 


k=\ 


-  l^akrxx(  -*)  +  °l  for    m  =  0 

for    m  <  0 


(2.18) 


k=\ 


The  relationship  between  the  autocorrelation  sequence  and  a  pure  moving  average 
model  may  be  written  by  setting  p  =  0  in  Eq.  (2.17) 


rxx(™)  = 


Je  /    bkb"\k-m) 

k=m 


for     m>  q 

for     0  <  m  < 
for     m  <  0 


(2.19) 


In  this  case  we  should  note  that 


K  =  bk       for    1  <  k  <  q 
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III.      TECHNIQUES  TO  SOLVE  FOR  THE  PARAMETERS  OF  ARMA 

MODELS 

The  estimation  of  the  parameters  of  ARMA  processes  is  a  classical  problem  which 
is  still  being  investigated  by  statisticians.  Methods  to  find  the  parameters  for  purely  AR 
processes  are  well  known,  but  for  ARMA  processes  some  problems  remain. 

There  is  a  nonlinear  relationship  between  the  ARMA  parameters  and  the 
autocorrelation  of  process  x(n).  The  nonlinear  equation  (2.17)  presents  the  difficulty  of 
estimating  the  ARMA  parameters,  even  when  we  know  the  autocorrelation  sequence 
exactly.  Techniques  based  on  iterative  maximum  likelihood  estimation  (MLE)  can  be 
used  to  find  the  ARMA  parameters.  These  techniques  require  complex  computations 
and  are  not  guaranteed  to  converge,  or  they  may  converge  to  the  wrong  solution. 
Therefore  they  are  not  practical  for  real  time  series.  For  AR  parameters,  techniques 
based  on  the  least  squares  criterion  lead  to  solutions  of  linear  equations  and  hence  re- 
duce the  computational  complexity.  Unfortunately,  the  moving  average  parameters  of 
an  ARMA  model  cannot  be  found  easily  by  solving  a  set  of  linear  equations.  The  MA 
parameters  are  convolved  with  the  impulse  response  coefficients  h(k)  which  causes  a 
nonlinear  relationship  between  the  autocorrelation  sequence  and  the  filter  coefficients. 

In  section  III-A  and  III-B,  we  will  discuss  the  methods  of  AR  and  MA  parameter 
estimation,  and  in  section  III-C  we  will  present  an  iterative  approach  to  find  the  pa- 
rameters. 

A.      AR  PARAMETER  ESTIMATION 

In  this  section,  we  present  three  widely  used  methods  of  extracting  the  model  pa- 
rameters from  a  given  block  of  measured  data  x(n). 
These  methods  are: 

1.  The  autocorrelation,  or  Yule-Walker  method 

2.  The  covariance  method 

3.  Burg's  method 

All  three  methods  of  estimating  AR  parameters  are  based  on  least-squares  minimization 
criteria  obtained  by  replacing  the  ensemble  averages  by  appropriate  time  averages  [Ref. 
7]- 
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The  criteria  for  the  optimal  forward  (  e-(n)  )  and  backward  (  e~(n)  )  predictors  are 
obtained  bv  minimizing 


£[»)2]         and 


E[e;{nf] 


where  e^{n)  and  e~(n)  are  the  result  of  filtering  x„  through  the  prediction-error  filter  as 
given  by 

eUn)  =  xn  +  apxn_}  +  apxn_2  + +  a  _* 


"pp~n-p 


ep  {n)  =  xn_p  +  apxn_p^  +apx„_p+2  + +  appxn 

The  autocorrelation  method  is  the  most  obvious  and  straightforward  one.  Equation 
(2.18)  may  be  evaluated  for  the  p+  1  lag  indices  0  <  m<p  and  put  into  the  following 
matrix  form 


(0) 


'XX 

'xx(l) 

'«(2) 


xx(-l) 


rA  -2) 

'-,,(-1) 
rxx(0) 


rxx(P)      rXx(P~l)      rJj-2) 


rxx(~P) 

1 

2 

o 

rxx(  -P  +  0 

a\ 

e 

0 

rxx(  -P  +  2) 

a2 

= 

0 

'xx(O) 

aP 

0 

(3.1) 


Bv  noting  that 


rxx(  ~k)  =  rxx(k) 


Eq.  (3.1)  can  be  written  as 


'xx(D  ^(0) 

r«(2)      fjji) 


rxx(2) 
rAO) 


rxxip)      r^ip-l)      r^ip-2) 


rxx(P) 
rxx(P~  1) 

1 

a, 

2 

0 

^O7  -  2) 

«2 

= 

0 

'xx(O) 

Op 

0 

(3.2) 


From  Eq.  (3.2) 


°\  =  ZJjakfa      where 


an  =  1 


i=0 
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We    can    replace    the    autocorrelations    rxx(k)    by    the    corresponding    sample 
autocorrelations  (biased  estimate)  computed  from  the  given  block  of  data,  where 

A-l-fc 

for      0<k<p  (3.3) 


We  use  the  biased  autocorrelation  estimator  because  the  unbiased  autocorrelation 
estimates  may  result  in  autocorrelation  matrices  that  are  not  positive  semidefmite,  which 
means  certain  matrix  equations  have  no  solution.  On  the  other  hand,  the  autocorrelation 
matrices  formed  from  the  biased  autocorrelation  estimate  will  always  be  positive  semi- 
defmite [Ref.   2  ]. 

By  solving  the  normal  equations,  we  can  obtain  estimates  of  the  model's  parameters 
{au  a2,  a3,. ...,ap,  a2,}.  Equation  (3.1)  is  known  as  the  AR,  Yule-Walker  or  Normal 
Equations  .  The  autocorrelation  matrix  of  this  equation  is  both  Toeplitz  and  Hermitian 
because 

?xxi  ~k)  =  r*xx(k) 


The  solution  of  the  Hermitian  Toeplitz  equations  can  be  computed  with  the 
Levinson  Algorithm.    This  algorithm  is  a  lattice  realization  for  linear  prediction  filters. 
We  illustrate  the  prediction-error  sequence  in  Figure  4. 

The  pth  prediction  error  is  given  by 

ep{n)  =  xn  +  apxn_x  +  apxn_2  + +  appx„_p  (3.4) 

Looking  at  the  first  two  prediction  errors 

el(n)  =  xn  +  anxn_] 

e2{n)  =  xn  +  a21xn_^  +  a22xn_2 

where  (1,  an),  and  (1,  a2i,  a22)  represent  the  best  predictors  of  orders  p=  1  and  p=  2,  re- 
spectively.  The  extra  index  is  used  to  indicate  the  order  of  the  predictor. 

In  the  autocorrelation  method  the  ensemble  average  in  minimization  of  the  forward 
prediction  error  is  replaced  by  the  least-squares  time  average  criteria  q>  as  follows 
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Figure  4.       Prediction-error  Sequence 


A'+r-i 

A  1        V     I    +/   \|2 

,V+/>-1 


H=0 


-r/i  +  /  ,aplcxn-i 


k=\ 


(3.5) 


where  it  is  assumed  that  the  data  xQ,  „y,, :%_,  are  observed  and  jc„  =  0  for  outside  the 

range  0  <  /»  <  Ar  —  1.  The  minimization  of  the  time-average  criteria  with  respect  to  the 
real  and  imaginary  part  of  the  ark  's  will  lead  exactly  to  the  same  set  of  Yule-Walker 
Equations  (3.1).  One  way  to  solve  these  equations  is  via  the  Levinson  recursion  which 
is  an  iterative  technique  that  extracts  the  next  order  predictor  from  the  previous  one. 
The  Levinson  algorithm  can  be  summarized  as  follows  [Ref.  7): 
1-  Initialize  the  recursion  at  p  =  0,  by  setting 


A0(z)  =  1 


and         Eq  =  rxx{0)  =  £[^] 


where  A0{z)  is  the  prediction-error  filter  and  £„  is  the  mean-squared  predic- 
tion error  at  the  zero  stage.   So,  initially  we  have  no  prediction.  At  stage  p, 
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If 


the  prediction-error  filter  will  be  Ap(z)  and  the  mean-squared  prediction  er- 
ror will  be  Ep. 

Compute  Yp+1  which  is  called  the   reflection   or  PARCOR  coefficient: 

p 

aplr{p  +1-/) 

/=o 

v.= — -P (3-6) 

ZvO 

/=0 

3-  Recursively  determine  the  (p  +  l)r/!  order  prediction-error  filter  polynomial 
Ap+X{z)\ 

Ap+](z)  =  Ap(z)  -  Y/)+1z-<p+1)^(z-1)  (3.7) 

4-  Update  the  mean-squared  prediction  error: 

Ep+l  =  (l-Y2p+l)Ep  (3.8) 

5-  Continue  the  iteration  until  the  final  desired  order  is  reached. 

If  the  process  x(n)  is  AR(p),  then  iteration  will  continue  up  to  order  p.  It  will  pro- 
vide the  AR  coefficients  alp,  a2p app  which  are  also  the  best  prediction  coefficients. 

If  we  continue  iteration  after  p,  all  prediction  coefficients  of  order  higher  than  p  will  be 
close  to  zero  [Ref.  7]. 

Although  the  autocorrelation  method  is  the  most  obvious  and  efficient  one,  and  the 
resulting  prediction-error  filter  is  guaranteed  to  be  minimal  phase,  it  suffers  from  the 
effect  of  prewindowing  the  data  sequence  x(n)  by  padding  it  with  zeros  to  the  left  and 
to  the  right.  This  reduces  the  accuracy  of  the  method,  especially  when  we  have  short 
data  records. 

In  the  Covariance  method  the  ensemble  average  in  minimization  of  the  forward  pre- 
diction error  is  replaced  by  the  time  average  as  follows 


A'-l 


^~*  k=\ 


(3.9) 


n=p 

The  only  difference  between  this  method  and  autocorrelation  method  is  in  the  limits 
of  the  summation.  In  the  covariance  method  we  observe  all  the  data  points  needed  to 
compute  (p.    Since  we  do  not  need  the  data  outside  the  range  0  <  n  <  X  —  1  we  do  not 
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need  to  set  xn  explicitely  to  zero.  For  that  reason  the  covariance  method  seems  to  be 
more  realistic. 

If  we  perform  the  minimization,  we  find  the  AR  parameter  estimates  as  the  solution 
of  equations  (3.10)  [Ref.  5]. 


Cjcx(U)        ^(1,2) 

c«(2,l)      ^(2,2) 
cxxiP^      cxx{p,2) 


CjcxO./O 

~a\~ 

"Cjoc(l,0)" 

cxx(2,P) 

a2 

— 

Cx*(2,0) 

Cxx(P>P) 

_aP_ 

cxx(Pfi) 

(3.10) 


where 


(V-l 


c(/.*)  = 


N 


I 


xn-jxn-k 


n=p 


Note  that  cxx(J,k)  is  an  estimate  of  rxx(J  —  k)  .  The  c„(j,k)  uses  the  sum  of  only  N-p 
lag  products  to  estimate  the  autocorrelation  function  for  each  lag  even  more  lags  are 
available.  In  contrast  in  the  estimation  of  rxx(0)  the  autocorrelation  method  uses  all  data 
point,  while  the  covariance  method  uses  only  N-p  data  points  in  the  summation.  The 
minimization  gives  us  a  non-Toeplitz  cxx{j,k)  matrix.  This  implies  that  wre  can  not  use 
the  Levinson  algorithm  to  solve  Eq.  (3.10).  The  equations  may  be  solved  by  using  the 
Cholesky  decomposition  which  will  be  computationally  more  expensive.  The  estimated 
poles  using  this  method  are  not  guaranteed  to  lie  within  the  unit  circle.  As  an  example, 
consider  a  first  order  predictor  (p=  1)  and  a  length-three  (N=  3)  sequence  [Ref.  7].  The 
time-averase  criteria  will  be 


<P 


=  -j  /^  I  et(")  1 2  =  y  (C*i  +  a\  i*o)2  +  (*2  +  a\  \x\?) 


n=.\ 


When  we  minimize  q>  by  differentiating  it  with  respect  to  an  and  setting  the  deriv- 
ative to  zero  ,  we  have 

fo  +  an*0)*0  +  (*2  +  a\\x\)x\  =  ° 


«n=- 


x\x0  +  X2X\ 

2    .       2 
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Note  the  denominator  does  not  depend  on  the  variable  x2  and  if  x2  is  large  enough 
we  might  have  magnitude  of  au  greater  than  one.  In  this  case  we  do  not  get  a 
minimum-phase  prediction-error  filter  which  is  sometimes  not  desirable  in  practice. 

Burg's  method  estimates  the  reflection  coefficients  and  then  uses  the  Levinson 
recursion  to  estimate  the  AR  parameters.  Burg's  minimization  criteria  is  to  minimize 
both  the  forward  prediction  error  ep(n)  and  the  backward  prediction  error  e~(n): 


A'-l 


<p  = 


1 


X 


£[>;(,7)2  +  ,;(«)2] 


n=p 


The  computational  steps  are  summarized  below  [Ref.  7]: 
1-  Initialize,  by  setting  ; 

eo"(«)  =  e<T(w)  =  xn  f°r         0<n<  X -  I 

and 


(3.11) 


£o  = 


,V-1 


n=0 


At  stage  p-1,  the  prediction-error  filter  will  be  /4,_,(z)  which  is  the  Z  trans- 
form of  the  sequence  {1,  ap_u,  ap_ia, ,ap_Up_x)  .    The  mean-squared  error 

will  be  £,_,.  e-_x{n),      and     e~_x{n)  can  be  calculated  for  p  —  I  <  n  <  X  —  1. 


Compute  the  reflection  coefficient: 

A'-l 


0<?H(fl)fp_i(" 


-1) 


*,-, 


>?=/? 


(3.12) 


I< 


n=p 


etAn)2  +  ev_An-  1) 


p-p 


To  guarantee  the  minimum-phase  property,  Y,  must  have  a  magnitude  less 
than  one. 

Compute  Ap(z)  using  the  Levinson  recursion  given  by: 


1 

2p,\ 

2p,2 


AP,P-\ 
ap,p 


0 


7>-U 

3/>-l,2 


« 


■p-1^-1 


-Y, 


a 


a 


P-i-P-i 
0 


p-\,p-2 


*p-],\ 


(3.13) 
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4-  Compute  e^(n)  and  e~(n)  for  p  <n<  N  —  1 . 

ep(n)  =  V-i(«  -  1)  ~  V^-it") 

5-  Update  the  mean-squared  error  as  follows: 

E^d-Y/E^  (3.15) 

6-  Continue  the  iteration  until  p  equals  the  model  order. 

The  Burg's  method  estimates  the  poles  which  are  on  or  inside  the  unit  circle.  This 
is  due  to  the  property  |  YJ  <  1.  Therefore,  care  must  be  taken  to  deal  with  the  situation 
when  |  YJ  =  1,  as  this  causes  the  prediction  filter  to  become  non-minimum  phase. 

B.     MA  PARAMETER  ESTIMATION. 

The  most  obvious  approach  to  estimate  the  MA  parameters  would  be  to  solve  the 
nonlinear  equation  (2.19)  using  the  autocorrelation  sequence.  Solutions  of  Eq.  (2.19) 
involve  difficult  spectral  factorization  techniques  [Ref.  8]. 

There  is  another  approach  called  Durbin's  method  which  is  related  Kolmogorov's 
theorem  [Ref.  4]  and  is  based  on  a  high  order  AR  approximation  of  the  MA  process. 
The  AR  process  allows  results  using  only  linear  operations.  Let 

<? 
B{z)  =  1  +  Yj^'X 

represent  the  system  transfer  function  of  an  MA(q)  process,  and 


A00(z)=l+Yjakz~k 


k=\ 

represent  the  system  transfer  function  of  an  AR{oo)  process  that  is  equivalent  to  the 
MA(q)  process.  Therefore,  we  have 

B{z)AJz)=\  (3.16) 

The  inverse  Z-transform  of  the  Eq.  (3.16)  is: 


■J 

E.,    .       [1       for       m  =  0  ,_  ,_. 

W* -*(*)- JO       for       m=l,2,...,<7  <3-17) 


«=i 


where  a0  =  1  and  aA  =  0  for  /c  <  0.  Therefore,  the  MA  parameters  can  be  determined  from 
the  infinite- order  AR  model  by  solving  (3.17). 

In  practice,  one  can  calculate  high-order  AR(M)  parameters,  where  M>  q  .  Based 

on  these  parameter  estimates  (\,aM(\),aM(2), ,aM(M))  ,  an  error  in  the  MA  part  is 

computed  [Ref.  2]. 

<Wm)  =  aM{m)  +  2^bnaM{m  -  n)  (3.18) 

According  to  Eq.  (3.17)  the  error  should  be  zero  for  all  m  except  for  m=  0.  But  in 
practice,  the  error  will  not  be  zero  when  using  finite  data,  so  MA  parameter  estimates 
are  obtained  by  the  minimization  of  the  squared  error  power,  given  by 


-£*! 


Ml2 
M (3,19) 

This  estimation  procedure  is  an  approximate  maximum  likelihood  estimate  (MLE). 
The  approximate  MLE  procedure  using  Durbin's  method  for  MA  parameters  results  in 
the  following  estimation  [Ref.  5]. 

b  =  -karaa  (3-20) 

where 

M-  \l-J\ 

\_Raa]ij=  M+  !       y      anan+{Hl  for  i,j=  1,2, ,q 

71=0 

M-l 


\jaa~]i  =    A/+  j     /  fnK+i  for  /  —  1,2, 


n=0 
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Again  the  Levinson  algorithm  may  be  used  to  solve  Eq.  (3.20)  for  the  b  parameters. 
The  estimated  zeros  of  B(z)  will  be  inside  the  unit  circle  by  the  minimum-phase  property 
of  the  autocorrelation  method. 

In  summary,  Durbin's  method  first  uses  the  data  x0,  xu  ..., .%_,  to  find  a  large  order 
AR(M)  model  using  the  autocorrelation  method.  Then  using  these  AR  parameter  esti- 

A  A  A 

mates  (1,  a„  a2,...,  aM)  as  the  data,  (bit  b2, ...,  bq)  is  found. 

Another  technique  for  estimating  the  MA  parameters  is  to  use  the  MA  spectral 
estimator  [Ref.  5]  which  is  given  by  Eq.  (3.21)  based  on  sample  correlation  values. 


PMA(f)=Yj^k)e'J2*/]C  (3-2I) 

k=-q 

Since  theoretically  Eq.  (3.21)  is  equal  to  a]\  B{J)\2,  the  MA  parameters  can  be  found 
by  using  the  Spectral  factorization  theorem  [Ref.  9].  This  theorem  shows  that  any  ra- 
tional power  spectral  density  of  a  stationary  signal  xn  can  be  factored  into  a  minimum- 
phase  form 

Pxx(z)  =  c]B{z)B{z-')  (3.22) 

C.     ITERATIVE  ARMA  PARAMETER  ESTIMATION 

Let  us  consider  the  modelling  of  ARMA  transfer  function  as  given  in  Figure  5. 
The  computational  steps  are  summarized  below  for  the  iteration  assuming  know- 
ledge of  the  model  order: 

1-  Form  the  sample  correlation  of  the  time  series  y(n)  by  using  the  biased 
autocorrelation  estimator. 

2-  Get  the  AR  coefficients  using  the  correlation  lags  >  q  (to  minimize  the  MA 
influence). 

3-  Inverse  filter  the  original  time  series  y(n)  via  the  inverse  of  the  AR  filter  (use 
the  AR  coefficients  of  step  2)  to  get  x2(n) 

4-  Form  the  sample  correlation  of  the  time  series  x2{n)  by  using  the  biased 
autocorrelation  estimator. 

5-  Get  the  MA  coefficients  using  the  sample  correlation  of  x2{n)  . 

6-  Inverse  filter  the  original  time  series  y(n)  via  the  inverse  of  the  MA  filter 
(use  MA  coefficients  which  are  estimated  in  step  5)  to  get  Jc,(«).  We  assume 
that  the  MA  coefficients  are  of  minimum  phase  (all  roots  are  inside  the  unit 
circle). 
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■ 

Figure  5.       Modelling  of  ARMA  transfer  function. 
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Form  the   sample  correlation   of  time   series  xx{n)   by  using   the   biased 
autocorrelation  estimator. 

Get  the  AR  coefficients  using  all  lags  or  lags  >  q  depending  on  the  method 
used. 

Inverse  filter  the  original  y(n)  series  using  the  AR  coefficient  estimates 
which  are  obtained  in  step  8,  to  get  x2(n). 

Get  the  sample  correlation  of  x2{n)  by  using  the  biased  autocorrelation  es- 
timator. 

Compute  the  MA  coefficients  using  the  sample  correlation  of  x2(n)  . 

Compute  the  error  in  the  AR  and  the  MA  parts  as  given  below 

P  <7 

error{j-\)=\(Ji-a}~')2+\(b]i-bJ~')2     for  j  =  2,3,...,10(3.23) 


i=i 


/=] 


where  superscript  j  is  used  to  denote  the  number  of  iteration.  The  upper 
count  of  j  is  experimentally  chosen.  Iteration  continue  up  to  10  if  the  co- 
efficients do  not  converge. 

If  error  >  ).,  then  go  to  step  6,  else  exit  the  program  using  the  last  updates 
of  the  filter  coefficients.  If  j>  10  terminate  with  an  error  message.  Note  ). 
is  a  small  experimentally  chosen  number. 
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The  AR  coefficients  can  be  obtained  via  a  pseudoinverse  from  the  modified  Yule- 
Walker  equations  [Ref.  2].  Equation  (2.17)  may  be  rewritten  for  the  p  lag  indices 
q+  1  <m<q  +  p  ,  and  put  into  matrix  form 


rxx(<1+  1) 


'xx(<7) 


rAq  +  P+l)      rJ.q+p  +  2) 


rxx(q-p+  i) 

~«r 

rxx(<1  -  P  +  2) 

«2 

= _ 

>xx(?) 

A 

'xx(<7  +  1) 

rxx(<7  +  2) 

^x(<?  +  p) 


(3.24) 


We  can  compute  the  autocorrelation  sequence  for  lags  q-p+  1  to  q  +  p,  therefore 
the  AR  parameters  are  found  as  the  solution  of  Eq.  (3.24)  where  the  MA  parts  influence 
is  minimal.  The  autocorrelation  matrix  is  of  Toeplitz  form. 

The  Moore-Penrose  pseudoinverse  A=  of  an  m  x  n  matrix  A  provides  the  minimum- 
norm  least  squares  solution   to  Ax  =  b.   The  solution  is  given  by 

x  =  A*b 

where  x  is  a  n  x  1  vector  that  simultaneously  minimizes  the  squared  equation  error,  for 
a  given  m  x  1  vector  b  [Ref.  2].  If  m=  n  and  rank  A  is  n  (i.e.,  A  is  nonsingular),  then  the 
pseudoinverse  becomes  the  square  matrix  inverse  A-  =  A~l.  If  m>n  (i.e.,  more 
equations  than  unknowns)  and  rank  A  is  n,  then 


A*  =  (AHAr]AH 


x  =  {AHA)~xAHb 


and 


The  superscript  H  is  used  to  denote  the  Hermitian  transpose  operation.  This  is  the 
least  squares  solution  for  a  set  of  overdetermined  equations. 

In  step  2  of  the  iteration,  we  will  use  correlation  lags  greater  than  q  while  in  step 
8  either  all  lags  or  lags  greater  than  q  are  used  depending  on  the  method. 

For  VIA  coefficients  estimation  the  Cholesky  decomposition  is  used.  If  the  matrix 
A  is  square  and  Hermitian,  then  the  usual  triangular  factorization  takes  on  the  special 
form 


A  =  RR 


u 


(3.25) 
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where  R  is  a  lower  triangular  matrix  with  nonzero  real  principal  diagonal  elements.  This 
decomposition  is  called  the   Cholesky  decomposition  [Ref.  2]. 

For  the  VIA  part,  the  statistical  autocorrelation  is  given  by 

q-  \k\ 

*«(*)  =  °l  E  W«  (3.26) 

Let  B  be  the  lower  triangular  Toeplitz  matrix  where  the  matrix  elements  bni  are  given 
in  terms  of  the  impulse  response  of  the  filter  B(z)  [Ref.  7]: 

bM  =  bn_i 
and  let  the  autocorrelation  matrix  of  x„  be 

Rxx(i,J)  =  RxxV  -J) 
Then,  the  transpose  matrix  BH  will  have  matrix  elements 

(B%  =  bi_n 

and  Eq.  (3.26)  can  be  written  as 

q-  \H\ 
Rxx('lJ)  =  Rxx(l  -J)  =°\     \    h+l-Pn 


n=0 


Thus,  in  matrix  notation 


R^iiJ)  =  al(BBH)y 


Rxx  =  o]BBH  (3.27) 


which  is  related  to  Eq.  (3.25)  with  the  assumption  that 


o\=\ 


Therefore,  an  approximation  of  the  MA  parameters  can  be  found  from  the  Cholesky 
decomposition  of  the  correlation  matrix  Rxx. 
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IV.      SIMULATION  RESULTS 

In  this  chapter,  computer  simulation  results  are  presented  to  show  how  three  dif- 
ferent estimation  methods  work  on  various  ARMA  models.  Two  ARMA(2,2)  models 
which  differ  in  pole-zero  locations,  an  ARMA(2,3)  model  and  an  ARMA(3,4)  model  are 
used  as  test  models.  In  addition  two  realizations  for  each  ARMA  process  are  utilized. 
Two  hundred  data  points  are  used  in  computing  the  sample  autocorrelation  values.  The 
three  different  AR  estimation  methods  are  explained  below. 

Method  1 

The  AR  coefficients  are  obtained  via  a  square  matrix  inverse  using  the  modified 
Yule-Walker  equations.  Correlation  lags  greater  than  q  are  used  to  minimize  the  MA 
part  influence  at  the  first  calculation.  For  the  remainder  of  the  iteration  the  same  cor- 
relation lags  are  used  to  minimize  potential  influence  from  the  MA  part.  A  Cholesky 
factorization  is  used  to  find  the  MA  part  coefficients. 

Method  2 

The  AR  coefficients  are  obtained  via  a  square  matrix  inverse  using  the  modified 
Yule-Walker  equations.  Correlation  lags  greater  than  q  are  used  to  minimize  the  MA 
part  influence  at  first  calculation.  For  the  remainder  of  the  iteration  the  correlation  lags 
starting  from  zero  are  used  assuming  the  MA  part  contribution  has  effectively  been  re- 
moved. The  Cholesky  factorization  is  used  to  find  the  MA  part  coefficients. 

Method  3 

The  AR  coefficients  are  obtained  via  a  pseudoinverse  instead  of  a  square  matrix 
inverse  using  the  correlation  function  starting  at  the  zero  lag  after  first  calculation.  This 
allows  the  use  of  all  important  correlation  lags.  The  Cholesky  factorization  is  used  to 
find  the  MA  part  coefficients. 

In  addition,  observation  noise  is  added  to  one  of  the  ARMA  models  and  the  re- 
sulting noisy  sequence  is  processed  via  the  three  different  methods.  The  observation 
noise  is  independent  of  the  driving  noise.  It  is  white  Gaussian  noise  with  zero  mean  and 
unit  variance.   The  signal  to  noise  ratio  (SNR)  is  about  15  dB. 

The  computer  program  computes  the  differential  errors  in  the  AR  and  MA  parts 
as  given  in  Eq.(3.23),  and  then  compares  it  with  a  small  experimentally  established  value 
X  (i.e.,  X  =  0.0001).    If  the  error  is  less  than  X  the  program  is  terminated.  If  the  error  is 
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larger  than  X  the  program  is  reentered  at  step  6.  Also,  if  the  coefficients  do  not  converge 
the  program  will  stop  after  nine  iterations. 

Comparisons  of  the  true  and  estimated  coefficient  differences,  of  pole-zero  lo- 
cations, of  distances  between  the  true  and  estimated  pole-zero  locations  and  of  radial 
differences  between  the  true  and  estimated  pole-zero  locations  are  presented  to  show 
how  well  the  three  methods  work.  Also  the  spectra  of  the  ARM  A  models  are  plotted  by 
using  the  true  and  estimated  coefficients. 

A.      THE  ARMA(2,2)  MODEL-A 

The  pole-zero  locations  for  this  model  are  illustrated  in  Figure  6. 


Figure  6.       The  ARMA(2.2)  Model-A,  Pole-zero  Locations 

1.      METHOD  1 

a.      Noise  Realization  I 

The  coefficients  converge  after  two  iterations.  Tables  1  and  2  present  the 
results. 
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Table  1.  ARMA(2,2)  MODEL-A,  METHOD  1, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

a, 

-1.80 

-1.8412 

-0.0412 

a2 

0.85 

0.8972 

+  0.0472 

K 

1.00 

1.3053 

+  0.3053 

bx 

0.80 

0.6439 

-0.1561 

b2 

0.80 

0.0572 

-0.742S 

Table  2.     ARMA(2,2)  MODEL-A,  METHOD  1,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.9  +  0.2j 

0.9206  +  0.2230J 

0.0308 

0.0189 

Pi 

0.9  -  0.2j 

0.9206  -  0.2230] 

0.030S 

0.0189 

-1 

-0.4  +  0.8j 

-0.3771 

0.8003 

1.1071 

Z2 

-0.4  -  0.8] 

-0.1162 

0.8488 

1.1071 

b.      Noise  Realization  2 

Using  a  different  noise  realization,  the  coefficients  still  converge  after  two 
iterations.  Tables  3  and  4  present  the  results. 


Table  3.  ARMA(2,2)  MODEL-A,  METHOD  1, 
COEFFICIENTS  COMPARISON, 
REALIZATION  2 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.80 

-1.8712 

-0.0712 

a2 

0.85 

0.9167 

+  0.0667 

h 

1.00 

1.1083 

+  0.1083 

bx 

0.80 

0.5071 

-0.2929 

h 

0.80 

0.1502 

-0.6498 
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Table  4.     ARMA(2,2)  MODEL-A,  METHOD  1,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.9  +  0.2] 

0.9356  +  0.2033] 

0.0358 

0.0047 

Pi 

0.9  -  0.2j 

0.9356  -  0.2033] 

0.0358 

0.0047 

Z\ 

-0.4  +  0.8j 

-0.2288  +  0.2884J 

0.5394 

0.2070 

Z2 

-0.4  -  0.8j 

-0.2288  -  0.2884J 

0.5394 

0.2070 

The  Spectra  using  the  true  and  the  estimated  network  coefficients  are 
plotted  in  Figure  7. 
2.      METHOD  2 

a.      Noise  Realization  1 

The  coefficients  converge  after  two  iterations.    Tables  5  and  6  present  the 
results. 


Table  5.  ARMA(2,2)  MODEL-A,  METHOD  2, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.80 

-2.0065 

-0.2065 

a2 

0.85 

1.0544 

+  0.2044 

K 

1.00 

1.3866 

+  0.3866 

ft, 

o.so 

0.6467 

-0.1533 

b. 

0.80 

-0.0178 

-0.8178 
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Table  6.     ARMA(2,2)  MODEL-A,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.9  +  0.2j 

1.0032  +  0.2189] 

0.1049 

0.0038 

P2 

0.9  -  0.2j 

1.0032-0.2189J 

0.1049 

0.0038 

Zl 

-0.4  +  0.8j 

-0.4924 

0.8053 

1.1071 

h 

-0.4  -  0.8j 

0.0261 

0.9064 

2.0344 

b.      Noise  Realization  2 

Using  a  different  noise  realization,  the  coefficients  still  converge  after  two 
iterations.  Tables  7  and  8  present  the  results. 


Table  7.  ARMA(2,2)  MODEL-A,  METHOD  2, 
COEFFICIENTS  COMPARISON, 
REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.S0 

-1.9267 

-0.1267 

a2 

0.85 

0.9729 

+  0.1229 

h 

1.00 

1.0982 

+  0.09S2 

bx 

0.80 

0.4598 

-0.3402 

b2 

0.80 

0.0902 

-0.7098 

Table  8. 


ARMA(2,2)  MODEL-A,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON. REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.9  +  0.2j 

0.9634  +  0.21 1  Sj 

0.0644 

0.0022 

P2 

0.9  -  0.2j 

0.9634- 0.21  lSj 

0.0644 

0.0022 

Zl 

-0.4  +  0.8j 

-0.2094  +  0.1958J 

0.6335 

0.3553 

Z2 

-0.4  -  O.Sj 

-0.2094-  0.1958] 

0.6335 

0.3553 
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The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  8. 

3.      METHOD  3 

a.      Noise  Realization  1 

The  coefficients  converge  after  two  iterations.   Tables  9  and  10  present  the 
results. 


Table  9.  ARMA(2,2)  MODEL-A,  METHOD  3, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

0i 

-1.80 

-1.8375 

-0.0375 

a2 

0.85 

0.S94S 

+  0.0448 

b0 

1.00 

1.3062 

+  0.3062 

by 

o.so 

0.6482 

-0.1518 

k 

0.80 

0.0657 

-0.7343 

Table   10.     ARMA(2,2)    MODEL-A,    METHOD    3,    POLE-ZERO    COM- 
PARISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

A 

0.9  +  0.2j 

0.9 1S7  +  0.2253] 

0.0314 

0.0218 

ft 

0.9  -  0.2j 

0.9187  -  0.2253J 

0.0314 

0.0218 

2i 

-0.4  +  O.Sj 

-0.3542 

0.8013 

1.1071 

h 

-0.4  -  0.8j 

-0.1421 

0.8405 

1.1071 

b.      Noise  Realization  2 

For  a  different  noise  realization,  the  coefficients  still  converge  after  one  it- 
eration. Tables  11  and  12  present  the  results. 
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Table  11.  ARMA(2.2)  MODEL-A,  METHOD 
3,  COEFFICIENTS  COMPAR- 
ISON, REALIZATION  2 


Coeff. 

True 

Estimated 

Difference 

<*\ 

-1.80 

-1.8545 

-0.0545 

<h 

0.85 

0.9043 

+  0.0543 

b0 

1.00 

1.1144 

+  0.1144 

b, 

0.80 

0.5272 

-0.2728 

b2 

o.so 

0.1  SOS 

-0.6192 

Table  12. 

ARMA(2,2)    MODEL-A,    METHOD 
PAR1SON,  REALIZATION  2 

3,    POLE-ZERO    COM- 

Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

A 

0.9  +  0.2j 

0.9272  +  0.21 1  lj 

0.0293 

0.0051 

Pi 

0.9  -  0.2j 

0.9272  -  0.21 1  lj 

0.0293 

0.0051 

zi 

-0.4  +  O.Sj 

-0.2365  +  0.3260J 

0.5014 

0.1639 

22 

-0.4  -  O.Sj 

-0.2365  -  0.3260J 

0.5014 

0.1639 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  9. 

In  summary,  the  estimated  coefficients  converge  after  two  iterations  in  all 
three  methods.  For  both  noise  realizations,  the  third  method  gives  the  best  result  in 
terms  of  the  coefficients  and  pole-zero  locations.  For  the  noise  realization  2,  the  coeffi- 
cients converge  after  one  iteration  using  method  3. 

In  each  method  the  AR  part  coefficients  of  the  ARMA(2,2)  model  tend  to 
be  more  accurate  than  the  MA  part  coefficients. 

B.      THE  ARMA(2,2)  MODEL-B 

The  pole-zero  locations  for  this  model  are  illustrated  in  Figure  10. 
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Figure  7.       The  Spectra  of  ARMA(2,2)  Model-A,  Method  1 
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AKMA(2,2)  Model  A,  Method  2.  Realization   1 
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Figure  8.       The  Spectra  of  ARMA(2.2)  Model-A,  Method  2 
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Figure  9.       The  Spectra  of  ARMA(2,2)  Model-A,  Method  3 
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The  ARMA(2,2)  Model-B,  Pole-zero  Locations 


1.      METHOD  1 

a.      Noise  Realization  1 

The  coefficients  converge  after  three  iterations.  Tables  13  and  14  present 
the  results. 


Table  13.  ARMA(2,2)  MODEL-B,  METHOD 
1,  COEFFICIENTS  COMPAR- 
ISON, REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

a, 

-1.20 

-1.2614 

-0.0614 

a2 

0.52 

0.4143 

-0.1057 

b0 

1.00 

1.0518 

+  0.0518 

bx 

0.50 

0.4045 

-0.0955 

b2 

0.3125 

-0.1317 

-0.4442 
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Table   14.     ARMA(2,2)  MODEL-B,  METHOD  1,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.6  +  0.4j 

0.6307  4-  0.1285J 

0.2732 

0.3870 

Pi 

0.6  -  0.4j 

0.6307-0.1285] 

0.2732 

0.3870 

Z] 

-0.25  +  0.5j 

-0.5950 

0.6074 

1.1071 

?2 

-0.25  -  0.5j 

0.2104 

0.6796 

2.0344 

b.      Noise  Realization  2 

For  a  different  noise  realization,  the  coefficients  converge  after  two  iter- 
ations. Tables  15  and  16  present  the  results. 


Table  15.  ARMA(2,2)  MODEL-B,  METHOD 
1,  COEFFICIENTS  COMPAR- 
ISON. REALIZATION  2 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.20 

-1.6704 

-0.4704 

a2 

0.52 

0.7879 

+  0.2679 

K 

1.00 

1.1018 

+  0.1018 

b> 

0.50 

0.1188 

-0.3812 

b2 

0.3125 

-0.350S 

-0.6633 

Table   16. 


ARMA(2,2)  MODEL-B,  METHOD  1,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.6  +  0.4j 

0.8352  +  0.3006] 

0.2553 

0.2425 

Pi 

0.6  -  0.4j 

0.8352  -  0.3006J 

0.2553 

0.2425 

zi 

-0.25  +  0.5j 

-0.6207 

0.6224 

1.1071 

~2 

-0.25  -  0.5j 

0.5129 

0.9121 

2.0344 
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The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  1 1. 

2.      METHOD  2 

a.      Noise  Realization  1 

The  coefficients  converge  after  six  iterations.   Tables  17  and  18  present  the 
results. 


Table  17.  ARMA(2,2)  MODEL-B,  METHOD 
2,  COEFFICIENTS  COMPAR- 
ISON. REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.20 

-1.1531 

+  0.0469 

a: 

0.52 

0.5013 

-0.0 1S7 

b0 

1.00 

1.0174 

+  0.0174 

b, 

0.50 

0.4732 

-0.0268 

b2 

0.3125 

0.1202 

-0.1923 

Table  18.     ARMA(2,2)  MODEL-B,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON. REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.6  +  0.4j 

0.5766  +  0.41  lOj 

0.0258 

0.0312 

P: 

0.6  -  0.4j 

0.5766-0.4110] 

0.0258 

0.0312 

zi 

-0.25  +  0.5j 

-0.2326  +  0.2532] 

0.2474 

0.2793 

z2 

-0.25  -  0.5j 

-0.2326  -  0.2532] 

0.2474 

0.2793 

b.      Noise  Realization  2 

For  a  different  noise  realization,  the  coefficients  converge  after  five  iter- 
ations. Tables  19  and  20  present  the  results. 
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Table  19.  ARMA(2.2)  MODEL-B,  METHOD 
2,  COEFFICIENTS  COMPAR- 
ISON. REALIZATION  2 


Coeff. 

True 

Estimated 

Difference 

a, 

-1.20 

-1.5914 

-0.3914 

<h 

0.52 

0.7540 

+  0.2340 

K 

1.00 

1.0610 

+  0.0610 

bx 

0.50 

0.1301 

-0.3699 

b2 

0.3125 

-0.2894 

-0.6019 

Table  20.     ARMA(2,2)  MODEL-B,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON. REALIZATION  2 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.6  +  0.4j 

0.7957  +  0.3477] 

0.2025 

0.1760 

P: 

0.6  -  0.4j 

0.7957  -  0.3477J 

0.2025 

0.1760 

*M 

-0.25  +  0.5j 

-0.5S72 

0.6030 

1.1071 

Zj 

-0.25  -  0.5j 

0.4645 

0.8720 

2.0344 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  12. 

3.      METHOD  3 

a.      Noise  Realization  1 

The  coefficients  converge  after  four  iterations.    Tables  21  and  22  present 
the  results. 
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Table  21.  ARMA(2.2)  MODEL-B,  METHOD 
3,  COEFFICIENTS  COMPAR- 
ISON, REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.20 

-1.1818 

+  0.0182 

«2 

0.52 

0.5245 

+  0.0045 

bo 

1.00 

1.0030 

+  0.0030 

*, 

0.50 

0.4436 

-0.0564 

h 

0.3125 

0.1015 

-0.2110 

Table  22.     ARMA(2,2)  MODEL-B,  METHOD  3,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

A 

0.6  +  0.4] 

0.5909  +  0.4187J 

0.0207 

0.0284 

Pi 

0.6  -  0.4j 

0.5909  -0.4187] 

0.0207 

0.0284 

~\ 

-0.25  +  0.5] 

-0.2211  +  0.2286] 

0.2729 

0.3050 

^2 

-0.25  -  0.5j 

-0.2211  -0.2286] 

0.2729 

0.3050 

b.      Noise  Realization  2 

For  a  different  noise  realization,  the  coefficients  converge  after  two  iter- 
ations. Tables  23  and  24  present  the  results. 


Table  23.  ARMA(2,2)  MODEL-B,  METHOD 
3,  COEFFICIENTS  COMPAR- 
ISON, REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

a, 

-1.20 

-1.5105 

-0.3105 

a2 

0.52 

0.6962 

+  0.1762 

K 

1.00 

1.048S 

+  0.04S8 

b 

0.50 

0.1827 

-0.3173 

h 

0.3125 

-0.2154 

-0.5279 
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Table  24.     ARMA(2.2)  MODEL-B,  METHOD  3,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

A 

0.6  +  0.4j 

0.7553  +  0.3547J 

0.1617 

0.1489 

Pi 

0.6  -  0.4j 

0.7553  -  0.3547J 

0.1617 

0.1489 

*i 

-0.25  +  0.5j 

-0.5486 

0.5823 

1.1071 

Z-2 

-0.25  -  0.5j 

0.3744 

0.7999 

2.0344 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  13. 

In  summary,  under  each  noise  realization  the  method  3  gives  the  best  result 
in  terms  of  the  coefficients  and  pole-zero  locations. 

For  each  method,  the  spectrum  of  the  ARMA(2,2)  model  using  estimated 
coefficients  indicates  reasonably  accurate  pole  locations  but  the  zero  locations  of  the 
spectrum  do  not  follow  the  original  ones.  This  means  that  the  AR  part  coefficients  of 
the  ARMA  model  tend  to  be  more  accurate  than  the  MA  part  coefficients. 

C.      THE  ARMA(3,4)  MODEL 

The  pole-zero  locations  for  this  model  are  illustrated  in  Figure  14. 
1.      METHOD  1 

a.      Noise  Realization  1 

The  coefficients  do  not  converge  in  nine  iterations.  Because  of  oscillations 
about  two  values,  the  average  of  the  two  values  is  used  as  an  estimate  of  the  coefficients 
in  the  comparison.   Tables  25  and  26  present  the  results. 
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Figure  1 1.       The  Spectra  of  ARMA(2,2)  Model-B,  Method  1 
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Figure  13.       The  Spectra  of  ARMA(2,2)  Model-B,  Method  3 
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Table  25.  ARMA(3,4)  MODEL,  METHOD  I, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coelf. 

True 

Estimated 

Difference 

<h 

-1.60 

-1.5275 

+  0.0725 

<h 

2.18 

1.5642 

-0.6158 

<h 

-1.36 

-1.0165 

+  0.3435 

(h 

0.7225 

0.4656 

-0.2569 

fro 

1.00 

3.3059 

+  2.3059 

*>, 

0.40 

1.3473 

0.9473 

i>2 

0.48 

-1.6482 

-2.1282 

l-h 

-0.32 

-2.4187 

-2.0987 
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Table  26.     ARMA(3,4)  MODEL,  METHOD   1,  POLE-ZERO  COMPAR- 
ISON. REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.2  +  0.9j 

0.0718  +  0.8101J 

0.1565 

0.1302 

Pi 

0.2  -  0.9j 

0.0718  -0.S101J 

0.1565 

0.1302 

Pi 

0.6  +  0.7j 

0.6919  +  0.4745J 

0.2435 

0.2610 

A 

0.6  -  0.7j 

0.6919  -  0.4745] 

0.2435 

0.2610 

a\ 

-0.4  +  O.Sj 

-0.6754  +  0.5652] 

0.3619 

0.4103 

Zi 

-0.4  -  0.8] 

-0.6754  -  0.5652] 

0.3619 

0.4103 

z3 

0.4 

0.9433 

0.5433 

0.0000 

b.      Noise  Realization  2 

Using  a  different  noise  realization,  the  coefficients  still  do  not  converge  in 
nine  iterations.  Because  of  oscillations  about  two  values,  the  average  of  the  two  values 
is  used  as  an  estimate  of  the  coefficients  in  the  comparison.  Tables  27  and  28  present 
the  results. 


Table  27.  ARMA(3,4)  MODEL,  METHOD  1, 
COEFFICIENTS  COMPARISON, 
REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.60 

-2.0479 

-0.4479 

a2 

2.18 

2.2690 

+  0.0890 

a3 

-1.36 

-1.4743 

-0.1143 

a4 

0.7225 

0.946S 

+  0.2243 

b0 

1.00 

2.6176 

+  1.6176 

h, 

0.40 

0.8737 

+  0.4737 

b2 

0.48 

-1.2400 

-1.7200 

bz 

-0.32 

-1.7033 

-1.3833 
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Table  28.     ARMA(3.4)  MODEL,  METHOD   1,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.2  +  0.9j 

0.0518  +  0.8258J 

0.1657 

0.1560 

Pi 

0.2  -  0.9j 

0.0518  -0.8258J 

0.1657 

0.1560 

Pz 

0.6  +  0.7j 

0.9722  +  0.6618] 

0.3741 

0.2644 

P* 

0.6  -  0.7j 

0.9722  -  0.6618] 

0.3741 

0.2644 

*\ 

-0.4  +  0.8j 

-0.6316  +  0.5489] 

0.3415 

0.3916 

z2 

-0.4  -  0.8j 

-0.6316  -  0.5489] 

0.3415 

0.3916 

h 

0.4 

0.9294 

0.5294 

0.0000 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  15. 

2.      METHOD  2 

a.      Noise  Realization  1 

The  coefficients  converge  after  nine  iterations.  Tables  29  and  30  present  the 
results. 


Table  29.     ARMA(3,4)  MODEL,  METHOD  2, 
COEFFICIENTS 
COMPARISON.REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

a, 

-1.60 

-1.6186 

-0.0186 

a2 

2.18 

2.2093 

+  0.0293 

«3 

-1.36 

-1.3948 

-0.0348 

a4 

0.7225 

0.7366 

+  0.0141 

K 

1.00 

1.1463 

+  0.1463 

bx 

0.40 

0.2664 

-0.1336 

b2 

0.48 

-0.0262 

-0.5062 

h 

-0.32 

-0.3386 

-0.0186 
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Table  30.     ARMA(3,4)  MODEL,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.2  +  0.9j 

0.2033  +  0.9077J 

0.0083 

0.0016 

Pi 

0.2  -  0.9j 

0.2033  -  0.9077] 

0.0083 

0.0016 

P> 

0.6  +  0.7j 

0.6060  +  0.6957] 

0.0073 

0.0079 

Pa 

0.6  -  0.7j 

0.6060  -  0.6957J 

0.0073 

0.0079 

z\ 

-0.4  +  0.8j 

-0.4197  +  0.5572J 

0.2435 

0.1819 

h 

-0.4  -  0.8j 

-0.4197-0.5572] 

0.2435 

0.1819 

^3 

0.4 

0.6070 

0.2070 

0.0000 

b.      Noise  Realization  2 

Using  the  second  noise  realization,  the  coefficients  converge  after  nine  it- 
erations. Tables  31  and  32  present  the  results. 


Table  31.     ARMA(3.4)  MODEL,  METHOD  2, 

COEFFICIENTS      COMPARISON, 

REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.60 

-1.6668 

-0.0668 

a2 

2.18 

2.3645 

+  0.1S45 

«3 

-1.36 

-1.5282 

-0.1 6S2 

a4 

0.7225 

0.S527 

+  0.1302 

b0 

1.00 

0.9923 

-0.0077 

b} 

0.40 

0.1274 

-0.2726 

b2 

0.48 

0.2222 

-0.2578 

b> 

-0.32 

-0.1221 

+  0.1979 
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Table  32.     ARMA(3.4)  MODEL,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.2  +  0.9j 

0.2063  +  0.9315J 

0.0321 

0.0007 

P2 

0.2  -  0.9j 

0.2063  -0.9315] 

0.0321 

0.0007 

Pi 

0.6  +  0.7j 

0.6271  +  0.7372] 

0.0460 

0.0037 

A 

0.6  -  0.7j 

0.6271  -  0.7372] 

0.0460 

0.0037 

£] 

-0.4  +  0.8] 

-0.22S6  +  0.5674J 

0.2889 

0.0806 

h 

-0.4  -  0.8j 

-0.2286  -  0.5674J 

0.2889 

0.0806 

Z3 

0.4 

0.32S7 

0.0713 

0.0000 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  16. 

3.      METHOD  3 

a.      Noise  Realization  1 

The  coefficients  converge  after  six  iterations.  Tables  33  and  34  present  the 
results. 


Table  33.  ARMA(3,4)  MODEL,  METHOD  3, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

a, 

-1.60 

-1.5S21 

-0.0179 

a2 

2.18 

2.1352 

-0.0448 

«3 

-1.36 

-1.3329 

+  0.0271 

aA 

0.7225 

0.7023 

-0.0202 

b. 

1.00 

1.1733 

+  0.1733 

b> 

0.40 

0.3140 

-0.0880 

b2 

0.48 

-0.0604 

-0.5404 

b~, 

-0.32 

-0.3978 

-0.0778 
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Table  34.     ARMA(3,4)  MODEL,  METHOD  3,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P. 

0.2  +  0.9j 

0.1913  +  0.9j 

0.0087 

0.0092 

Pi 

0.2  -  0.9] 

0.1913  -0.9j 

0.0087 

0.0092 

Pi 

0.6  +  0.7j 

0.5997  +  0.6855] 

0.0145 

0.0101 

Pi 

0.6  -  0.7] 

0.5997  -  0.6855] 

0.0145 

0.0101 

zi 

-0.4  +  0.8j 

-0.4533  +  0.5691] 

0.2369 

0.20S9 

2; 

-0.4  -  0.8] 

-0.4533  -  0.5691] 

0.2369 

0.20S9 

•^3 

0.4 

0.6406 

0.2406 

0.0000 

b.      Noise  Realization  2 

Using  the  second  noise  realization,  the  coefficients  converge  after  four  it- 
erations. Tables  35  and  36  present  the  results. 


Table  35.  ARMA(3,4)  MODEL,  METHOD  3, 
COEFFICIENTS  COMPARISON, 
REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

«i 

-1.60 

-1.5671 

+  0.0329 

a2 

2. IS 

2.2062 

+  0.0262 

«3 

-1.36 

-1.4005 

-0.0405 

a* 

0.7225 

0.7728 

+  0.0503 

K 

1 .00 

0.9932 

-0.0068 

b} 

0.40 

0.2049 

-0.1951 

b2 

0.48 

0.2119 

-0.2681 

b. 

-0.32 

-0.1654 

+  0.1546 
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Table  36.     ARMA(3,4)  MODEL,  METHOD  3,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.2  +  0.9j 

0.1876  +  0.9336J 

0.0358 

0.0203 

Pi 

0.2  -  0.9j 

0.1876  -  0.9336J 

0.0358 

0.0203 

Pi 

0.6  +  0.7j 

0.5959  +  0.7051] 

0.0065 

0.0069 

P* 

0.6  -  0.7j 

0.5959  -0.7051] 

0.0065 

0.0069 

h 

-0.4  +  0.8j 

-0.2937  +  0.5923] 

0.2333 

0.0033 

%2 

-0.4  -  0.8j 

-0.2937  -  0.5923] 

0.2333 

0.0033 

^3 

0.4 

0.3810 

0.0190 

0.0000 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plot- 
ted in  Figure  17. 

When  using  method  1  the  coefficients  do  not  converge.  They  oscillate 
about  two  sets  of  values  in  an  alternating  fashion,  hence  the  average  of  the  two  sets  is 
used  as  the  estimate  of  the  coefficients  in  the  comparison.  The  parameters  of  the  first 
realization  converge  after  nine  iterations  using  method  2,  and  after  six  iterations  using 
method  3.  The  parameters  of  the  second  realization  converge  after  four  iterations  using 
method  3.  Method  3  provides  the  best  results  in  terms  of  the  coefficients  and  pole-zero 
locations  for  both  noise  realizations.  The  coefficients  converge  also  earlier  when  using 
method  3  compared  to  the  other  two  methods.  Method  2  performs  not  as  well  but  gives 
also  reasonable  results. 

The  spectrum  due  to  the  poles  of  the  ARMA(3,4)  models  using  the  esti- 
mated coefficients  closely  resembles  the  original  spectrum.  The  AR  part  coefficient  es- 
timates are  more  accurate  than  the  MA  part  coefficient  estimates. 

D.      THE  ARMA(3,4)  MODEL  WITH  OBSERVATION  NOISE 

Observation  noise  is  added  to  the  ARMA(3,4)  model  and  resulting  noisy  sequence 
is  processed  via  the  three  different  methods.  The  observation  noise  is  independent  of  the 
driving  noise.  It  is  white  Gaussian  noise  with  zero  mean  and  unit  variance.  The  signal 
to  noise  ratio  is  about  15  dB. 

This  model  is  illustrated  in  Figure  18. 
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Figure  15.       The  Spectra  of  ARMA(3,4)  Models,  Method  1 
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ARMA(3,4)  Model,  Method  2.  Realization   1 
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Figure  16.       The  Spectra  of  ARMA(3,4)  Models,  Method  2 
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ARMA(3,4)  Model,  Method  3,  Realization   1 
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Figure  17.       The  Spectra  of  ARMA(3,4)  Models,  Method  3 


52 


Observation     Noise 
w(n) 


Driving  Noise 
£(n) 


ARMA(3,4) 


x(n) 


y(n) 


Figure   18.       The  ARM A(3,4)  Model  with  Observation  Noise 

1.      METHOD  1 

The  coefficients  do  not  converge  in  nine  iterations.  Because  of  oscillations 
about  two  values,  the  average  of  the  two  values  is  used  as  an  estimate  of  the  coefficients 
in  the  comparison.   Tables  37  and  38  present  the  results. 


Table  37.  ARMA(3,4)  MODEL  WITH  OB- 
SERVATION NOISE,  METHOD  1, 
COEFFICIENTS  COMPARISON. 


Coeff. 

True 

Estimated 

Difference 

«] 

-1.60 

-2.0285 

-0.4285 

a3 

2.18 

1.9241 

-0.2559 

a3 

-1.36 

-1.4285 

-0.0685 

(h 

0.7225 

0.8145 

+  0.0920 

b„ 

1.00 

4.5732 

+  3.5732 

bi 

0.40 

1.7315 

+  1.3315 

b% 

0.48 

-2.4570 

-2.9370 

b3 

-0.32 

-3.3648 

-3.0448 
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Table  38.     ARMA(3,4)      MODEL      WITH      OBSERVATION      NOISE, 
METHOD  1,  POLE-ZERO  COMPARISON 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.2  +  0.9j 

0.0147  +  0.8350] 

0.1963 

0.2010 

Pi 

0.2  -  0.9] 

0.0147-0.8350J 

0.1963 

0.2010 

Pi 

0.6  4-  0.7j 

0.9995  +  0.4109] 

0.4931 

0.4721 

P* 

0.6  -  0.7j 

0.9995  -  0.4109] 

0.4931 

0.4721 

zi 

-0.4  +  0.8j 

-0.6723  +  0.5565] 

0.3652 

0.4157 

z2 

-0.4  -  0.8] 

-0.6723  -  0.5565] 

0.3652 

0.4157 

~3 

0.4 

0.9660 

0.5660 

0.0000 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plotted 
in  Figure  19. 

2.      METHOD  2 

The  coefficients  do  not  converge  in  nine  iterations.  Because  of  oscillations 
about  two  values,  the  average  of  the  two  values  is  used  as  an  estimate  of  the  coefficients 
in  the  comparison.   Tables  39  and  40  present  the  results. 


Table  39.  ARMA(3,4)  MODEL  WITH  OB- 
SERVATION NOISE,  METHOD  2, 
COEFFICIENTS  COMPARISON. 


Coeff. 

True 

Estimated 

Difference 

«] 

-1.60 

-2.6027 

-1.0027 

a2 

2.1S 

3.8865 

+  1.7065 

a3 

-1.36 

-2.8006 

-1.4406 

<k 

0.7225 

1.4933 

+  0.7708 

K 

1.00 

S.5983 

+  7.5983 

bx 

0.40 

3.4234 

+  3.0234 

b2 

0.4S 

-3.2534 

-3.7334 

b> 

-0.32 

-4.6219 

-4.3019 
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Table  40.     ARMA(3,4)      MODEL      WITH      OBSERVATION      NOISE, 
METHOD  2.  POLE-ZERO  COMPARISON 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.2  +  0.9j 

0.3264  +  0.8586J 

0.1330 

0.1446 

Pi 

0.2  -  0.9j 

0.3264  -  0.8586] 

0.1330 

0.1446 

Pi 

0.6  +  0.7j 

0.9749  +  0.9052J 

0.4273 

0.1138 

P* 

0.6  -  0.7j 

0.9749  -  0.9052] 

0.4273 

0.1138 

zi 

-0.4  +  0.8j 

-0.6152  +  0.5170J 

0.3555 

0.4082 

z2 

-0.4  -  0.8] 

-0.6152  -  0.5170J 

0.3555 

0.4082 

z3 

0.4 

0.S323 

0.4323 

0.0000 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plotted 
in  Figure  20. 

3.      METHOD  3 

The  coefficients  converge  after  four  iterations.    Tables  41  and  42  present  the 
results. 


Table  41.  ARMA(3,4)  MODEL  WITH  OB- 
SERVATION NOISE,  METHOD  3, 
COEFFICIENTS  COMPARISON. 


Coeff. 

True 

Estimated 

Difference 

fli 

-1.60 

-1.4994 

+  0.1006 

flj 

2.1S 

2.1311 

-0.0489 

«3 

-1.36 

-1.3224 

+  0.0376 

<k 

0.7225 

0.7353 

+  0.0128 

K 

1.00 

1.3898 

+  0.3898 

*. 

0.40 

0.2875 

-0.1125 

b2 

0.48 

0.0298 

-0.4502 

b3 

-0.32 

-0.2462 

+  0.0738 
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Table  42.     ARMA(3,4)      MODEL      WITH      OBSERVATION      NOISE, 
METHOD  3,  POLE-ZERO  COMPARISON 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.2  +  0.9j 

0.1823  +  0.9331] 

0.0375 

0.0257 

Pi 

0.2  -  0.9j 

0.1823  -0.9331] 

0.0375 

0.0257 

Pi 

0.6  +  0.7] 

0.5674  +  0.7010] 

0.0326 

0.0281 

Pt 

0.6  -  0.7j 

0.5674  -  0.7010J 

0.0326 

0.0281 

zi 

-0.4  +  0.8j 

-0.3481  +  0.4908J 

0.3135 

0.1532 

?2 

-0.4  -  0.8j 

-0.3481  -0.4908] 

0.3135 

0.1532 

^3 

0.4 

0.4893 

O.OS93 

0.0000 

The  spectra  using  the  true  and  the  estimated  network  coefficients  are  plotted 
in  Figure  21. 

In  summary,  the  estimated  coefficients  do  not  converge  using  method  1  or  2 
but  they  do  converge  after  four  iterations  using  method  3.  The  spectra  using  method  1 
and  2  are  poor.  But  the  method  3  gives  results  close  to  the  actual  ones.  The  spectrum 
estimated  with  this  method  follows  the  original  pattern  except  for  the  zero  location. 
This  is  believed  to  be  partially  due  to  the  imprecise  MA  coefficient  estimation  procedure. 
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Figure  19.       The  Spectra  of  ARMA(3,4)  Models  with  Obser.  Noise,  Metliod  1 
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AKMA(3,4)  Model  with  Obser.  Noise,  Method  2 
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Figure  20.       The  Spectra  of  ARMA(3,4)  Models  with  Obser.  Noise,  Method  2 
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ARMA(3,4)  Model  with  Obser.  Noise.  Method  3 
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Figure  21.       The  Spectra  of  ARMA(3,4)  Models  with  Obser.  Noise,  Method  3 
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V.     CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  thesis,  an  iterative  technique  to  estimate  the  coefficients  of  ARMA  models 
driven  by  a  random  input  is  presented.  The  AR  parameters  are  estimated  initially,  then 
MA  parameters  are  estimated  assuming  the  AR  parameters  have  been  compensated  for. 
In  an  iterative  fashion  MA  and  AR  contributions  are  removed  from  the  original  data 
allowing  improved  AR  and  MA  coefficient  estimates.  Three  different  AR  estimation 
methods  are  experimentally  explored.  The  third  method  provides  the  best  result  in  terms 
of  the  true  and  estimated  coefficients  and  pole-zero  locations.  This  method  uses  the 
pseudoinverse  with  correlation  function  values  starting  at  the  zero  lag  after  removing  the 
VIA  influence.  For  models  of  real  time  data  which  have  an  odd  number  of  poles,  one 
should  address  the  issue  of  the  DC  component  (i.e.,  remove  it).  Simulation  results  for 
an  odd  ordered  pole  model  is  presented  in  Appendix  C. 

Also,  by  examining  the  spectra  of  the  models,  we  can  say  estimation  of  poles  is 
obtained  more  accurately  than  the  estimation  of  zeros.  The  Cholesky  factorization  is 
used  for  the  estimation  of  the  MA  part  coefficients  but  it  is  an  approximate  solution. 
For  that  reason,  we  do  not  expect  superior  results  for  the  VIA  estimation  part. 

Further  research  should  concentrate  on  improving  the  MA  part  coefficients  esti- 
mation. 
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APPENDIX  A.      COMPUTER  PROGRAM  TO  ESTIMATE 

COEFFICIENTS 

^'oV'sY'jV'jV'jV'sV's'fjVVr'iV  irk  3V}V}W«V3Wr?V}V?V}V}WoV?V}V'>Wc?V}V,5V,!V 

%  This  program  written  by  * 
%  Gurhan  Kayahan  for  IBM  PC/AT  * 
%     Mat lab  package.  * 

^*^VVr,sViViV,'5VyriV'!V<VVfiVVr,!VVc,Vc'iV'5V,!ViVVriViV?VVr*iV**iV*»VVrVf 
o/  **-  y-  »'-  -j u  -i-  j-  ~'-  »T-  j-  -'-  -J-  j-  -»-  -^  »'-  «*-  y -  -*-  y-  «*-  Jj  Jf  «'f  -•-  j-  y-  »*-  j-  ->-  »*-  j-  »'-  j u  Jf  j- 

%  GENERATION  OF  ARMA  TIME  SERIES   * 

0/  y-  y-  y-  -  ■  - y-  y- y-  y-  -  •  -  - u  y-  y-  y-  y-  y-  y-  - '  -  - '  -  y -  -  •  -  *  ■  - » '  -  y-  -  *  -  y-  y-  -  ■-  y-  y-  y-  -  *  -  y-  »u  y-  y* 

i=sqrt(-l); 

rootsl=C  -0.  25  +  0.  5i  -  0.  25  -  0.  5i];  change  zeros  for  each  model 

roots2=  [0.  6  +  0.  4i  0.  6  —  0. 4i  0. 8]; change  poles  for  each  model 

b=poly(rootsl); f ind  true  coefficients  of  MA  part 

a=poly(roots2); f ind  true  coefficients  of  AR  part 

a=real(a); 

t=l:  1000; 

y=t*0; 

rand( 'normal' ); generate  random  signal 

stddev=sqrt( 1); 

rnum=stdev'vrand(t); 

yl=f ilter(b, a, rnum);  filter  random  signal  through  ARMA  model 

x=  [ones{l  ,999)];  generate  step  signal 

y2=f ilter(b,a,x); find  step  response  of  ARMA  model  to  discard  transient  response 

y=yl(52:252); 

0/  y-  y-  y-  y-  y-  y-  JU  y-  y-  y-  y-  y~  y-  y-  y-  y-  j-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y  .  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  ju  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y-  y, 

%   CALCULATE  AR  COEFFICIENTS  USING  YULE -WALKER  APPROACH  * 

/o"  "  "  "  "  "  "  "  "  "  "  "  "  "  "  '*  "  "  "  "  "  "  "  "  '•  "  "  "  "  ?»  "  "  "  "  "  "  "  *«  "  "  "  "  "  "  'v  '•  •«  "  »  "  "  '•  "  "  «  «» 

rl=xcorr(y); 

rl=rl(201:  399); 

r=flipy(rl); 

row=r(194: 196); 

col=flipy(r(192: 194)); 

a=toeplitz( col , row); 

b=flipy(-r(191:  193)); 

c=b' ; 

x=b"inv(a);  square  matrix  inverse 

al=x'; 

%   CALCULATE  MA  COEFFICIENTS  USING  CH0LESKY  FACTORIZATION  * 

fo.s  ,*  ,*  ,„  *s  ,»  ,»  ,*  ,*  „  ,„  ,»  <s  ,%  „  ,,  4%  ,»  tx  ,%i*  ,»  „  «  „  ,»  "wwn  '*  *»  «**« 

d=  [1]; 

n=  [1  aid)  aU2)  <9l(3)]; 

z=f ilter(n,d,y); 

zl=z(4:201); 

rzl=xcorr(zl); 

rz=rzl(198:  395); 

i=l; 

while   i<199 

rz3(i)=rz(i)/198; 
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-  ~'-  ^f-  JU  ^-  JU  JU  J1.  J«  **-  .*-  J1* 


end 

rz2=rz3(  1:  3); 
rzt=toeplitz( rz2); 
cho=chol( rzt); 
bl=cho(l,:  ); 

%   ROUTINE  FOR  ESTIMATION  OF  AR  COEFFICIENTS  * 

i=2; 
while  i<ll 

d=  [^(i-1,1)  M(i-1,2)  M(i- 1,3)3; 

n=  [13; 

zar=filter(n,d,y); 

rar=xcorr(zar); 

rl=rar(201:  399); 

r=flipy(rl); 

row=r( 194: 196); change  correlation  lags  for  each  method 

col=flipy(r( 192: 194)); change  correlation  lags  for  each  method 

g=toeplitz(col ,row); 

e=flipy( -r( 191: 193) ); change  correlation  lags  for  each  model 

f=e'; 

x=f"'vinv(g);  square  matrix  inverse  (replace  with  x=pinv(g)  for  pseudo  inverse) 

a(i,:  )=x  ; 

al=  [al',  a(i,:  )']'; 

%   ROUTINE  FOR  ESTIMATION  OF  MA  COEFFICIENTS 

d=  [13; 

n=  [1  a(i,l)  ad, 2)   a(i,3)3; 

z=f ilter(n,d,y); 

zl=z(4:  201); 

rzl=xcorr(zl); 

rz=rzl(198:  395); 

k=l; 

while  k<199 

rz3(k)=rz(k)/198; 
k=k+l; 
end 

rz2=rz3(l:  3); 
rzt=toeplitz(rz2); 
cho=chol( rzt); 
bt(i,:  )=cho(l,:  ); 
bl=  Cbl',  bt(i,:  )'3'; 

e(i-l)=sum(  (al(i,:  )-al(i-l,:  )).  a  2  )/3+sum(  (bl(  i,:  )-bl(  i-1, :  ) ).  a  2  )/3; 
if  e(i-l)  <  0.0001 

al=al; 

bl=bl; 
else 

i=i+l; 
end 
end 
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APPENDIX  B.   COMPUTER  PROGRAM  TO  CALCULATE  TRUE  AND 

ESTIMATED  SPECTRA 

%   CALCULATION  OF  TRUE  SPECTRUM 

a=  [1  -1.20  0.52];  change  true  poles  for  each  models 

b=  [1  0.50  0.3125];  change  true  zeros  for  each  models 

al=  [<9  zerosil,  123)];  zero  padding 

bl=  lb  zeros(l,  124)];  zero  padding 

fal=fft(al);  fast  fourier  transform  of  AR  part 

fbl=fft(bl);  fast  fourier  transform  of  MA  part 

mfal=abs( fal).  a  2; 

mfbl=abs(fbl).  a  2; 

c=mfbl.  /mfal; 

a=max(c); 

i-i; 

while  i<129 

c(i)=c(i)/a; 

i=i+l; 
end 
spec=10*logl0(c); 

%   CALCULATION  OF  ESTIMATED  SPECTRUM 

a=  [1  -1.1531  0.5013];  change  estimated  poles  for  each  models 

b=  [1.0174  0.4732  0.1202];  change  estimated  zeros  for  each  models 

al=  la  zeros(  1,123)]; 

bl=  lb  zerosil,  124)]; 

fal=fft(al); 

fbl=fft(bl); 

mfal=abs( fal).  a  2; 

mfbl=abs(fbl).  a  2; 

c=mfbl.  /mfal; 

a=max(c); 

while  i<129 

c(i)=c(i)/a; 

i=i+l; 
end 
spec7=10*logl0(c); 

%  PLOTTING  OF  TRUE  AND  ESTIMATED  SPECTRUM 

O,'  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JU  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  -t-  JL  JL  JL  JL  -"-  JL  JL  JL  JL  JU  JU  JL  -'-  -'-  JL 
^**  *»  #»  **  *»  *t  n  *»  *\  /»  #*  *v  #»  rff  #»  **  ?f  #*  #\  rt  **  f*  **  *t  n  /»  *«  **  n  #»  *»  #»  /»  *»  #*  **  *»  »»  *»  *v  **  ?»  « 

t=0:  127; 

t=t'; 

plot(t(l:65),spec7(l:65),'*'  , t(  1:  65)  ,spec(  1:  65)) 

title('ARMA(2,2)  Model,  Method  3,  Realization  l') 

xlabel( ' Sample  Number') 

ylabel( 'Magnitude(dB) ' ) 

grid 

text(25,-5,' Using  true  coefficients') 
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text(25,-9, '****  After  4  iterations') 
delete  sped,  met 
meta  spec! 


64 


APPENDIX  C.      SIMULATION  RESULTS  FOR  AN  ARMA(2,3)  MODEL 

For  models  with  an  odd  number  of  poles,  one  of  the  poles  must  be  located  on  the 
real  axis  in  the  z-domain  to  generate  real  time  data. 

The  pole-zero  locations  for  this  model  are  illustrated  in  Figure  22. 
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Figure  22.       The  ARMA(2,3)  Model  Pole-zero  Locations 

1.      METHOD  1 

a.      Noise  Realization  1 

The  coelTicients  do  not  converge  in  nine  iterations.  Because  of  oscillations 
about  two  values,  the  average  of  the  two  values  is  used  as  an  estimate  of  the  coefficients 
in  the  comparison.    Tables  43  and  44  present  the  results. 
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Table  43.  ARMA(2,3)  MODEL,  METHOD  1, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

«i 

-2.00 

-1.9425 

+  0.0575 

a2 

1.48 

1.4028 

-0.0772 

«3 

-0.416 

-0.3639 

+  0.0521 

K 

1.00 

1.0562 

+  0.0562 

*, 

0.50 

0.5184 

-0.01S4 

b2 

0.3125 

0.1681 

-0.1444 

Table  44.     ARMA(2,3)  MODEL,  METHOD   1,  POLE-ZERO  COMPAR- 
ISON. REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

Pi 

0.6  +  0.4j 

0.6526  +  0.3809J 

0.0559 

0.0596 

Pi 

0.6  -  0.4j 

0.6526  -  0.3809J 

0.0559 

0.0596 

Pi 

0.8 

0.6373 

0.1627 

0.0000 

2\ 

-0.25  +  0.5j 

-0.2454  +  0.3145] 

0.1855 

0.1989 

z2 

-0.25  -  0.5j 

-0.2454  -  0.3145] 

0.1855 

0.1989 

b.      Noise  Realization  2 

Using  a  different  noise  realization,  the  coefficients  converge  after  seven  it- 
erations. Tables  45  and  46  present  the  results. 
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Table  45.  ARMA(2,3)  MODEL,  METHOD  1, 
COEFFICIENTS  COMPARISON, 
REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

a, 

-2.00 

-2.8571 

-0.8571 

<h 

1.48 

2.7525 

1.2725 

a3 

-0.416 

-0.9162 

+  1.3322 

K 

1.00 

1.3418 

+  0.3418 

bx 

0.50 

0.1291 

-0.3709 

b2 

0.3125 

-0.4401 

-0.7526 

Table  46.     ARMA(2.3)  MODEL,  METHOD   1,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

A 

0.6  +  0.4j 

0.S301  +  0.2765J 

0.2611 

0.2664 

Pi 

0.6  -  0.4j 

0.8301  -  0.2765] 

0.2611 

0.2664 

Pi 

0.8 

1.1968 

0.3968 

0.0000 

-M 

-0.25  +  0.5j 

-0.6229 

0.6237 

1.1071 

Zj 

-0.25  -  0.5j 

-0.5266 

0.9236 

2.0344 

2.      METHOD  2 

a.      Noise  Realization  1 

The  coefficients  do  not  converge  in  nine  iterations.  Because  of  oscillations 
about  two  values,  the  average  of  the  two  values  is  used  as  an  estimate  of  the  coefficients 
in  the  comparison.   Tables  47  and  48  present  the  results. 
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Table  47.  ARMA(2,3)  MODEL,  METHOD  2, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

«i 

-2.00 

-1.6534 

+  0.3466 

a2 

1.48 

0.9727 

-0.5073 

<h 

-0.416 

-0.2044 

+  0.2116 

K 

1.00 

1.2677 

+  0.2677 

*. 

0.50 

0.S274 

+  0.3274 

b2 

0.3125 

0.3741 

+  0.0616 

Table  48.     ARMA(2,3)  MODEL,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON. REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.6  +  0.4j 

0.5267  +  0.2515J 

0.1656 

0.1425 

Pi 

0.6  -  0.4j 

0.5267  -  0.2515] 

0.1656 

0.1425 

Pi 

0.8 

0.6001 

0.1999 

0.0000 

-i 

-0.25  +  0.5j 

-0.3263  +  0.4349] 

0.1006 

0.1807 

-: 

-0.25  -  0.5] 

-0.3263  -  0.4349] 

0.1006 

0. 1 807 

b.      Noise  Realization  2 

Using  a  second  noise  realization,  the  coefficients  do  not  converge.  Because 
of  oscillations  about  two  values,  the  average  of  the  two  values  is  used  as  an  estimate  of 
the  coefficients  in  the  comparison.   Tables  49  and  50  present  the  results. 
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Table  49.  ARMA(2,3)  MODEL,  METHOD  2, 
COEFFICIENTS  COMPARISON, 
REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

«i 

-2.00 

-1.1005 

+  0.8995 

a2 

1.48 

-0.1144 

-1.5944 

«3 

-0.416 

0.3945 

+  0.8105 

^0 

1.00 

1.8387 

+  0.8387 

bi 

0.50 

1.3064 

+  0.S064 

b2 

0.3125 

0.4440 

+  0.1315 

Table  50.     ARMA(2,3)  MODEL,  METHOD  2,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.6  +  0.4j 

0.8145  +  0.2883] 

0.2418 

0.2478 

Pi 

0.6  -  0.4j 

0.8145  -0.2883J 

0.2418 

0.2478 

fh 

0.8 

-0.5285 

1.3285 

3.1415 

~\ 

-0.25  +  0.5j 

-0.3553  +  0.3395] 

0.1919 

0.3444 

z2 

-0.25  -  0.5j 

-0.3553  -  0.3395] 

0.1919 

0.3444 

3. 


METHOD  3 

a.      Noise  Realization  1 

The  coefficients  converge  after  five  iterations.  Tables  51  and  52  present  the 


results. 
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Table  51.  ARMA(2,3)  MODEL,  METHOD  3, 
COEFFICIENTS  COMPARISON, 
REALIZATION  1. 


Coeff. 

True 

Estimated 

Difference 

«, 

-2.00 

-1.1057 

+  0.8943 

a2 

1.48 

0.2111 

-1.2689 

«3 

-0.416 

0.0890 

+  0.5050 

K 

1.00 

1.8594 

+  0.8594 

b, 

0.50 

1.4885 

-0.9885 

b2 

0.3125 

0.8320 

-0.5195 

Table  52.     ARMA(2,3)  MODEL,  METHOD  3,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  1. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.6  +  0.4j 

0.6493  +  0.2j 

0.2059 

0.2892 

P: 

0.6  -  0.4j 

0.6493  -  0.2j 

0.2059 

0.2S92 

Pi 

0.8 

-0.1928 

0.9928 

3.1415 

zi 

-0.25  +  0.5j 

-0.4002  +  0.5360] 

0.0509 

0.1777 

Z-> 

-0.25  -  0.5j 

-0.4002  -  0.5360] 

0.0509 

0.1777 

b.      Noise  Realization  2 

Using  a  second  noise  realization,  the  coefficients  do  not  converge.  Because 
of  oscillations  about  two  values,  the  average  of  the  two  values  is  used  as  an  estimate  of 
the  coefficients  in  the  comparison.   Tables  53  and  54  present  the  results. 
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Table  53.  ARMA(2,3)  MODEL,  METHOD  3, 
COEFFICIENTS  COMPARISON, 
REALIZATION  2. 


Coeff. 

True 

Estimated 

Difference 

«i 

-2.00 

-1.7973 

+  0.2027 

«2 

1.48 

0.9959 

-0.4841 

«3 

-0.416 

-0.9354 

-0.5194 

^0 

1.00 

1.2540 

+  0.2540 

6, 

0.50 

0.6150 

+  0.1150 

^ 

0.3125 

0.0151 

-0.2974 

Table  54.     ARMA(2,3)  MODEL,  METHOD  3,  POLE-ZERO  COMPAR- 
ISON, REALIZATION  2. 


Poles- 
Zeros 

True 

Estimated 

Distance 

Radial  Diff. 

P\ 

0.6  +  0.4j 

0.1263  +  0.7679] 

0.5997 

0.2564 

Pi 

0.6  -  0.4j 

0.1263  -0.7679J 

0.5997 

0.2564 

Pz 

0.8 

1.5446 

0.7446 

0.0000 

2\ 

-0.25  +  0.5j 

-0.4645 

0.5440 

1.1071 

Z2 

-0.25  -  0.5j 

-0.0259 

0.5479 

1.1071 

Simulation  results  have  shown  that  for  the  ARMA(2,3)  model  with  one 
pole  on  positive  real  axis,  realization  2  gives  poor  results.  This  is  because  the  DC  com- 
ponent for  this  realization  is  -  0.9017  while  it  is  -0.0580  for  realization  1.  For  parameter 
estimations  the  DC  component  should  be  removed.  Therefore  the  data  is  filtered  effec- 
tively removing  the  pole  on  the  positive  real  axis  (i.e.,  DC  component  is  removed  and 
the  ARMA(2,3)  model  becomes  ARMA(2,2)  model).  The  results  of  this  ARMA(2,2) 
model  was  presented  in  chapter  IV-B.  The  ARMA(2,3)  and  ARMA(2,2)  true  and  esti- 
mated spectra  using  method  1  are  presented  in  Figure  23.  Method  2  results  are  pre- 
sented in  Figure  24  and  method  3  results  are  presented  in  Figure  25.  The  simulation 
shows  that  the  method  3  gives  good  results  as  long  as  the  DC  component  is  small. 
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Figure  23.       The  Spectra  of  ARiMA(2,3)  and  ARMA(2,2)  Model-13,  Method  1 
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Figure  24.       The  Spectra  of  ARMA(2,3)  and  ARMA(2,2)  Model-B,  Method  2 
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